Phase behaviour and the random phase approximation 
for ultrasoft restricted primitive models 
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Phase separation of tiie ultrasoft restricted primitive model (URPM) with Gaussian charges is re- 
investigated in the random phase approximation (RPA) — the 'Level A' approximation discussed by 
Nikoubashman, Hansen and Kahl [J. Chem. Phys. 137, 094905 (2012)]. We find that the RPA pre- 
dicts a region of low temperature vapour-liquid coexistence, with a critical density much lower than 
that observed in either simulations or more refined approximations (we also remark that the RPA 
critical point for a related model with Bessel charges can be solved analytically). This observation 
suggests that the hierarchy of approximations introduced by Nikoubashman et al. should be analo- 
gous to those introduced by Fisher and Levin for the restricted primitive model [Phys. Rev. Lett. 71, 
3826 (1993)], which makes the inability of these approximations to capture the observed URPM 
phase behaviour even more worthy of investigation. 



Recently Coslovich, Hansen and Kahl (CHK) intro- 
duced a novel class of Gaussian charge cloud models for 
mixtures of interpenetrable polycations and polyanions 
in solution [1] [2] . The low temperature phase behaviour 
of these models was explored both by Monte-Carlo and 
molecular dynamics simulations ^ [2| , and in mean field 
theory by Nikoubashman, Hansen and Kahl (NHK) [3]. 
Our interest in this class of models stems from a different 
perspective. In mesoscale models, particularly in dissipa- 
tive particle dynamics (DPD) [4] soft interactions are the 
norm. Then it is both natural, and indeed essential, to 
smear out point charges into charge clouds. The diver- 
gence of the long-range Coulomb law as r — )■ (where r 
is the center-center separation) is replaced by a smooth 
cutoff, thus ensuring thermodynamic stability according 
to a theorem by Fisher and Ruelle [5] . The precise form 
of the charge smearing is often tuned to the numerical al- 
gorithm used to calculate the electrostatic interactions, 
and a consensus on the best approach has yet to emerge 
[51 [7] . Whilst for mesoscale modelling applications the 
low temperature phase behaviour is not in itself of pri- 
mary importance, the screening properties though are of 
great interest and our research into this aspect will be 
reported more thoroughly elsewhere. 

The canonical example of this class of models, which 
CHK termed the ultrasoft restricted primitive model 
(URPM), is an equimolar mixture of Gaussian charge 
clouds, which are identical apart from the sign of the 
charges, and for which only the electrostatic interactions 
are kept. The URPM is a natural counterpart to the well- 
studied restricted primitive model (RPM) of equi-sized 
charged hard spheres [SHTU] in which the short-range 
Coulombic divergence is hidden behind the hard core re- 
pulsion. For the URPM, CHK reported a region of low 
temperature vapour-liquid phase coexistence, for which 
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the terminus on increasing temperature bears many of 
the hallmarks of a tricritical point. Above this point, 
and a likely reason for the apparent tricriticality, is either 
a weak second order transition or a rapid crossover be- 
tween an insulating dielectric phase of neutral 'molecules' 
of paired opposite charges and a conducting 'plasma' 
phase containing a substantial fraction of free ions. Sub- 
sequently NHK investigated a hierarchy of mean-field ap- 
proximations in an attempt to understand in detail the 
origin of the low temperature URPM phase behaviour. 
This hierarchy was built in analogy to the earlier work of 
Fisher and Levin on the RPM [lljJJJ. The simplest level 
of approximation, termed 'Level A' in NHK, is analogous 
to the Fisher-Levin DH (Debye-Hiickel) approximation. 
It is identical to the random phase approximation (RPA) 
from integral equation theory, and incorporates the mu- 
tual attractions and repulsions in a linearised way. The 
next level of approximation, 'Level B' in Ref. [3] and 
DHBj (Debye-Hiickel-Bjerrum) in the Fisher-Levin clas- 
sification, captures the formation of ion pairs — a crucial 
aspect of the non-linear physics at low temperatures. 

NHK assert that "there is no phase separation at [the 
'Level A'] approximation" (below Eq. (30) in Ref. (2). 
This caught our attention, as we have known for some 
time that the RPA for a related Bessel charge model 
(discussed below) does exhibit phase separation, with 
a critical point which can be determined analytically. 
Prompted by this discrepancy, our own further investi- 
gations reveal that the RPA for Gaussian charges does 
have a region of phase separation, but at a much lower 
density than investigated by NHK. 

To set the problem up, we consider an equimolar mix- 
ture of — iV_ = N/2 charge clouds (polyions) in a 
volume V, with an overall density p = N/V. Gaussian 
charge clouds interact with the following pair potential, 



(1) 



where u{r) is the pair potential between charge clouds 
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of the same sign, /? = l/k-oT is the inverse of the tem- 
perature measured in units of Boltzmann's constant, 
is the Bjerrum length which plays the role of a coupling 
constant, r is the separation, and cr is a measure of the 
size of the charge cloud. For Gaussian charges the ra- 
dial charge distribution corresponding to this potential 
is (27rcr2)-3/2e-'''/2<T=. The function erf(r/2cr) - r as 
r — >■ 0, thus ensuring the Coulombic divergence is re- 
placed by a smooth cutoff. 

An interesting alternative to the Gaussian charge 
URPM is provided by a Bessel charge model. For this 
case the interaction potential is simply 



/3u(r) 



'(1 



(2) 



This corresponds to a radial charge distribution 
K\(r I a)l (^-K^u^r) where Ki is a modified Bessel func- 
tion (hence the name). Although this radial charge dis- 
tribution diverges as r — >■ 0, the interaction potential 
itself is again smoothly cutoff. 

Gaussian charges are blessed by being particularly well 
suited to the Ewald summation method for handling long 
range Coulomb interactions, as has been noted by CHK. 
Bessel charges, on the other hand, are not so well suited 
for simulations but provide perhaps the simplest non- 
trivial example of an ultrasoft primitive model when it 
comes to analytical work. Obviously, the definition of a 
in the two potentials cannot be exactly matched up and 
this should be born in mind when making comparisons. 

In reciprocal space these potentials are 



/3i2(fc) 

where, writing q = ka, 



w{q) = 



w{ka) 



(Gaussian), 



1/(1 + ^2) (Bessel). 



(3) 



(4) 



The definition of a in the two models is chosen to match 
up the long wavelength behaviour here. 

The random phase approximation (RPA) for this class 
of models takes the form c±± = —(3u±± for the direct 
correlation functions ^ [31 [T3HT5] . Because of the ab- 
sence of hard cores, the RPA is also equivalent to the 
mean spherical approximation (MSA). From the RPA, 
the total correlation functions, h±±{r) = ±/i(r), follow 
by inversion of the Ornstein-Zernike equations. In recip- 
rocal space the solution is 



hik) 



—4TrlBw{ka) 
fc2 -I- k'^w{ka) 



(5) 



In this fcp = 



^ 'inlBP is the square of the Debye wavevector. 
It follows from Eq. ^ that the density-density struc- 
ture factor is given by S]\[]sf{k) = 1 and, somewhat less 
trivially, the charge-charge structure factor is given by 
Szzik) = ky[k^ + k^w{ka)]. 
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FIG. 1. RPA free energy for Gaussian URPM at ct^t/Zb = 
0F/27 « 0.0656, from Eqs. with w{q) = e""'. A 

function Ap, with I3A = 4.03, is added to the free energy to 
reveal the common tangent construction without perturbing 
the phase behaviour. 



In all these we notice the prominent role played by the 
denominator D{k) = k'^ + k^w{ka). As is well known 
[21 E] the zeros of this function in the complex fc-plane 
determine the asymptotic behaviour of the total corre- 
lation functions, and are crucial to understanding the 
screening properties of the system particularly for ap- 
plications in mesoscale modelling. The asymptotic be- 
haviour typically crosses over from being purely expo- 
nential to being damped oscillatory as one increases the 
density past the so-called Kirkwood line in the density- 
temperature plane [TBj. More generally, this is referred 
to as a Fisher- Widom line [T7] . For Gaussian charges the 
asymptotic behaviour is determined by the complex roots 
of q^ + qY)e~'^ = 0, where go = ku<J. The most relevant 
roots are given by q^ = Wo{—q^) where Wq is the princi- 
pal branch of the Lambert W function [T5^. From this, or 
by direct calculations j3j , the Kirkwood line for Gaussian 
charges is given by qu = e"^/^ ~ 0.6065. The Kirkwood 
line for Bessel charges is determined by the complex roots 
of the biquadratic equation + + = 0. For qo < h 



the roots are all purely imaginary, whereas for go > 5 
they are all complex. Hence in this case the Kirkwood 
line takes the simple form go = |. 

Now we turn to the free energy. It follows from the 
density-density structure factor that the compressibility- 
route equation of state is trivially that of an ideal gas, 
for which the free energy density is 



/3f'i = p(lnip-l) 



(6) 



Note there are two species of ions contributing to this, 
each at a density ^p, and we have neglected the ther- 
mal de Broglie wavelength as it plays no role in phase 
coexistence. The virial-route equation of state, and the 
energy-route equation of state {via coupling constant in- 
tegration) give rise to the same result, which can be in- 
tegrated to a non-trivial excess free energy density. The 
result is 
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dq 
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{l + ^w{q))-qlwiq)]. (7) 
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FIG. 2. (color online) Vapour-liquid coexistence regions (bin- 
odals plus critical points) for the URPM with Gaussian or 
Bessel charges. Approximations are RPA (present work) and 
RPABj ('Level B' in Ref. 3 ). Data for RPABj is taken from 
Fig. 7 of Ref. |3] , and the simulation data is taken from Fig. 18 
of Ref. [2]. See Table [l] for locations of critical points. 



For point charges 'w{q) — 1 and this reduces to the exact 
DH limiting law /3/°'^ = ~k'^/12tt. The total free energy 
density, used in calculating the phase behaviour, is given 
by the sum of Eqs. ([6| and ([7| : 



/ = .r + r 



(8) 



In the Gaussian case Eq. ([7| is exactly equal to Eq. (29) 
in Ref. [3] . Figure [I] shows the total free energy, from 
Eq. ([s]) , as a function of density at a judiciously chosen 
temperature, illustrating the existence of a common tan- 
gent construction. The full phase behaviour is plotted in 
Fig. [2j marked 'RPA', where also are shown the 'Level 
B' results replotted from Ref. [J, here marked 'RPABj', 
and simulation results taken from Ref. [2]. The Gaus- 
sian RPA critical point, found numerically, is located at 
Ib/ct « 26.25 and pa^ « 1.014 x IQ-^ (see also Table 
|l]). This corresponds to a reduced Debye wavelength of 
(7d ~ 0.335 which places the critical point somewhat on 
the low density side of the RPA Kirkwood line. 

At this point we should comment on the choice of re- 
duced (dimensionless) temperature. CHK and NHK use 
Mo = u{0) as an energy scale but this frustrates direct 
comparison with the RPM. Our own preference is to use 
the long range behaviour of the potential characterised by 
the reduced Bjerrum length Ib/ct. Since (Huq — Ib/ct ^/tt 
for Gaussian charges, to facilitate the comparison with 
CHK and NHK we universally use /Ib as a reduced 
temperature. In this, cr is the parameter entering the in- 
teraction potentials in Eqs. ([T]) and ([2| for the URPM, 
and the hard sphere diameter for the RPM. 

For the Bessel case, the RPA excess free energy can 
be obtained in closed form. The last term in Eq. ([t]) 
evaluates to — ^□/(STrcr^) = — /Bp/(2cr). On multiplying 



system 


method 






Refs. 


RPM 


DH 


0.11 
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URPM (B) 


RPA 


0.085 
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TABLE L Vapour-hquid critical points for the RPM and 
URPM with Gaussian (G) and Bessel (B) charges. RPA re- 
sults are those reported in the present paper. For the Gaus- 
sian URPM the RPABj ('Level B' in Ref. [3]) and simulation 
results are taken from Table 1 in Ref. |^. All results (both 
here and in the main text) are accurate to the final digit. 



through by cr^, the first part of the integral is 



1 

4^ 



dqqHn[l+ 



(9) 



To solve this, we learn from the (Schwinger-)Feynman 
parameter trick [21 and rewrite it as 



1 

4^ 



dq I du 



q^l 



(10) 



Making for the time being the assumption that q^ < \ 
(so that we are on the low density side of the Kirkwood 
line), the g- integral can now be done, by the method of 
partial fractions, to get 



87r^7o 



du 



(11) 



where z = \/l — 4u (hence the temporary restriction on 
) • We note that du = —-^zdz, so the w-integral can also 
be done. After taking careful account of the integration 
limits, the final result for the free energy is 



2y2-(l + z)3/2_(l_z)3/2 l^pa^ 



24ttV2 



(12) 



where now z — -y/l — iq^ . Whilst this result has been 
derived for q^ < |, it holds by analytic continuation for 
all (7d . As one crosses the Kirkwood line from low to high 
density, z crosses over from being purely real to purely 
imaginary, so that 



V49d - 1 (9D > i) 



(13) 



Nevertheless the free energy remains purely real and is 
continuous across the Kirkwood line. (Note that the 
roots of g** -f q2 _|_ g2 _ Q given by q^ = —1 ± ^z.) 



4 



Like the Gaussian case, the RPA free energy for the 
Bessel case has a region of vapour-hquid phase coexis- 
tence at low densities and temperatures. The critical 
point can be found by solving ff^f /dp^ — d^f/dp^ — 
from Eqs. ([6|, ^ and (12). An analytic solution can 
be obtained, which is Ib/o^ = 12y/3 « 20.78 and pa'^ = 
l/(487r\/3) « 3.829 x lO'^ (see also Table |l]). This cor- 
responds to z = iVS and = 1, thus for Bessel charges 
the RPA critical point lies on the high density side of the 
Kirkwood line. The phase behaviour for the Bessel case, 
calculated numerically, is also shown in Fig. [2] 

Table |l] compares the vapour-liquid critical points for 
the RPM and the URPM, using various approximations. 
We see that the DH approximation for the RPM, and 
the RPA for the URPM, both predict critical points at 
low densities and temperatures. When Bjerrum pair- 
ing is incorporated {i. e. DHBj for RPM, and RPABj for 
Gaussian URPM), the critical temperature remains un- 
changed but the critical density is considerably increased. 
For the RPM, this brings the predicted critical point 
quite close to the simulations, within 20% for the critical 
temperature (for a detailed discussion, see Ref. [H]). For 
the URPM with Gaussian charges though, the predicted 
critical point is still considerably distant from the simu- 



lations. In particular the predicted critical temperature 
is at least a factor of three above the observed value. We 
can to some extent confirm this observation, as we have 
looked for phase separation in the Gaussian URPM using 
Monte-Carlo methods, at temperatures in the vicinity of 
the RPA critical point, and have found no evidence of 
such. This singular aspect of the phase behaviour of the 
URPM stands in marked contrast to the RPM. Some 
possible explanations have been proposed by NHK [3]. 

The observation that the critical temperature remains 
unchanged in comparing RPA and RPABj can be traced 
to the fact that in the latter approximation the Bjerrum 
pairs are an ideal spectator species piMT^ . As such they 
cannot, in themselves, influence the phase behaviour of 
the unpaired ions. The quasi-chemical equilibrium be- 
tween paired and unpaired ions changes the coexistence 
densities, in accordance with the law of mass action, but 
the critical temperature itself remains unaffected. 

To summarise, the Fisher-Levin hierarchy of approx- 
imations developed for the RPM can be pursued also 
for the URPM, with similar trends, indicating the two 
models should show similar phase behaviour. The fact 
that they do not deepens the mystery uncovered by 
Nikoubashman, Hansen and Kahl in Ref. [7 and clearly 
warrants further investigation. 
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